*
* $Id$
*
* $Log: fkdres.F,v $
* Revision 1.1.1.1  2002/06/16 15:18:36  hristov
* Separate distribution  of Geant3
*
* Revision 1.1.1.1  1999/05/18 15:55:15  fca
* AliRoot sources
*
* Revision 1.1.1.1  1995/10/24 10:20:05  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.45  by  S.Giani
*-- Author :
*=== dres =============================================================*
*                                                                      *
      SUBROUTINE FKDRES ( M2, M3, T1, U, EREC, LOPPAR, JFISS )
 
#include "geant321/dblprc.inc"
#include "geant321/dimpar.inc"
#include "geant321/iounit.inc"
*
*----------------------------------------------------------------------*
*                                                                      *
*  New version of DRES created  by A.Ferrari & P.Sala, INFN - Milan    *
*                                                                      *
*  Last change  on  10-apr-93  by  Alfredo Ferrari, INFN - Milan       *
*                                                                      *
*  Dres93: Dres91 plus the RAL fission model taken from LAHET thanks   *
*          to R.E.Prael                                                *
*  Dres91: new version from A. Ferrari and P. Sala, INFN - Milan       *
*          This routine has been adapted from the original one of the  *
*          Evap-5 module (KFA - Julich). Main modifications concern    *
*          with kinematics which is now fully relativistic and with    *
*          the treatment of few nucleons nuclei, which are now frag-   *
*          mented, even though in a very rough manner. Changes have    *
*          been made also to other routines of the Evap-5 package      *
*                                                                      *
*----------------------------------------------------------------------*
*
*----------------------------------------------------------------------*
*                                                                      *
*  Input variables:                                                    *
*     M2 = Mass number of the residual nucleus                         *
*     M3 = Atomic number of the residual nucleus                       *
*     T1 = Excitation energy of the residual nucleus before evaporation*
*     U  = Excitation energy of the residual nucleus after evaporation *
*     Erec = Recoil kinetic energy of the residual nucleus             *
*            The recoil direction is given by Coslbr (i)               *
*                                                                      *
*  Significant variables:                                              *
*     JA = Present mass number of the residual nucleus                 *
*     JZ = Present atomic number of the residual nucleus               *
*     Smom1 = Energy accumulators for the six types of evaporated      *
*             particles                                                *
*                                                                      *
*    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     *
*    !!!! Please note that the following variables concerning !!!!     *
*    !!!! with the present residual nucleus must be set before!!!!     *
*    !!!! entering DRES91: Ammres, Ptres                      !!!!     *
*    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     *
*                                                                      *
*----------------------------------------------------------------------*
*
C---------------------------------------------------------------------
C SUBNAME = DRES --- EVAPORATION
C           EVAPORATION DATA SHOUD BE READ ON INPUT STAGE
C---------------------------------------------------------------------
C*****A LA EVAP III(TWA,8-68)
#include "geant321/eva0.inc"
#include "geant321/eva1.inc"
#include "geant321/forcn.inc"
#include "geant321/higfis.inc"
#include "geant321/inpflg.inc"
#include "geant321/labcos.inc"
#include "geant321/nucdat.inc"
#include "geant321/resnuc.inc"
      DIMENSION ZMASS (6), Z2MASS(6), C(3), Q(0:6), FLKCOU(6), CCOUL(6),
     &          THRESH(6), SMALLA(6), R(6), S  (6), SOS   (6), STRUN(6),
     &          EYE1  (6), EYE0  (6), SMOM1    (6), BNMASS(6)
      DIMENSION CORRRR(6)
      REAL RNDM(2)
      LOGICAL LOPPAR, PENBAR, LFIRST
      PARAMETER (EXPMIN=-88,EXPMAX=88)
      SAVE ZMASS, Z2MASS, EMHN, EMNUM, UM, AMUMEV, AMEMEV, QBRBE8,
     &     BNMASS, IEVEVP, NBE8, NRNEEP, LFIRST
      DATA IEVEVP / 0 /
      DATA LFIRST / .TRUE. /
*
      IEVEVP = IEVEVP + 1
C-------------------------------------- 1.ST CALL INIT
      IF ( LFIRST ) THEN
         LFIRST = .FALSE.
         FKEY   = ZERZER
         NBE8   = 0
         NRNEEP = 0
         EXMASS(1) = 1.D+03 * ( AMNEUT - AMUAMU )
         EXMASS(2) = FKENER ( ONEONE, ONEONE )
         EXMASS(3) = FKENER ( TWOTWO, ONEONE )
         EXMASS(4) = FKENER ( THRTHR, ONEONE )
         EXMASS(5) = FKENER ( THRTHR, TWOTWO )
         EXMASS(6) = FKENER ( FOUFOU, TWOTWO )
         ZMASS(1) = 1.D+03 * AMUAMU + EXMASS (1)
         ZMASS(2) = 1.D+03 * AMUAMU + EXMASS (2)
         ZMASS(3) = 2.D+03 * AMUAMU + EXMASS (3)
         ZMASS(4) = 3.D+03 * AMUAMU + EXMASS (4)
         ZMASS(5) = 3.D+03 * AMUAMU + EXMASS (5)
         ZMASS(6) = 4.D+03 * AMUAMU + EXMASS (6)
         BNMASS (1) = 0.D+00
         BNMASS (2) = 0.D+00
         BNMASS (3) = ZMASS (1) + ZMASS (2) - ZMASS (3)
         BNMASS (4) = TWOTWO * ZMASS (1) + ZMASS (2) - ZMASS (4)
         BNMASS (5) = ZMASS (1) + TWOTWO * ZMASS (2) - ZMASS (5)
         BNMASS (6) = TWOTWO * ( ZMASS (1) + ZMASS (2) ) - ZMASS (6)
         DO 1234 KK = 1,6
            Z2MASS (KK) = ZMASS (KK) * ZMASS (KK)
1234     CONTINUE
         AMUMEV = 1.D+03 * AMUAMU
         AMEMEV = 1.D+03 * AMELEC
         QBRBE8 = FKENER ( EIGEIG, FOUFOU ) - TWOTWO * EXMASS (6)
         EMN = 1.D+03 * AMNEUT
         EMH = ZMASS (2)
         TMP16 = 16.D+00
         UM  = AMUMEV + FKENER ( TMP16, EIGEIG ) / 16.D+00
         EMHN  = EMH - EMN
         EMNUM = EMN - UM
      END IF
*  |  End of initialization:
*  +-------------------------------------------------------------------*
C     --------------------------------- START OF PROCESS
*  +-------------------------------------------------------------------*
*  |  Initialize Npart and Smom if nothing has been already evaporated
*  |  for this event
      IF ( JFISS .LE. 0 ) THEN
         DO 775 I=1,6
            NPART(I) = 0
            SMOM1(I) = ZERZER
  775    CONTINUE
      END IF
*  |
*  +-------------------------------------------------------------------*
      JA = M2
      JZ = M3
      U  = T1
      RNMASS = 1.D+03 * AMMRES + U
* P2res and  Ptres are the squared momentum and the momentum of the
* residual nucleus (now in relativistic kinematics), Umo the
* invariant mass of the system!
      UMO  = RNMASS
      UMO2 = UMO * UMO
      ELBTOT  = RNMASS + EREC
      GAMCM   = ELBTOT / RNMASS
      ETACM   = 1.D+03 * PTRES / RNMASS
      HEVSUM  = ZERZER
 1000 CONTINUE
      LOPPAR = .FALSE.
*  +-------------------------------------------------------------------*
*  |                Check for starting data inconsistencies
      IF (JA-JZ .LT. 0) THEN
         WRITE(LUNOUT,6401)
         WRITE(LUNERR,6401)
 6401    FORMAT('1 Dres: cascade residual nucleus has mass no. less',
     &       ' than Z!!')
         RETURN
*  |
*  +-------------------------------------------------------------------*
*  |                Rough treatment for very few nucleon residual
*  |                nuclei. The basic ideas are:
*  |        a) as many as possible alpha particles are emitted
*  |        b) particles are emitted one per time leaving a residual
*  |           excitation energy proportional to number of nucleons
*  |           left in the residual nucleus (so we deal only with
*  |           two body kinematics)
*  |       T A K E   I N T O   A C C O U N T   T H A T   T H I S
*  |       T R E A T M E N T   I S   E X T R E M E L Y   R O U G H
*  |       T H E   T A S K   B E I N G   O N L Y   T O   S U P P L Y
*  |       S O M E T H I N G   T O   S H A R E   E N E R G Y   A N D
*  |       M O M E N T U M   A M O N G   A   F E W   F R A G M E N T S
      ELSE IF ( JA .LE. 6 .OR. JZ .LE. 2 ) THEN
*  | 1000 continue moved above according to FCA suggestion
*1000    CONTINUE
         JRESID = 0
         IF ( JA .GT. 4 ) GO TO 2000
*  |  +----------------------------------------------------------------*
*  |  | First check we are not concerning with a couple of neutrons or
*  |  | protons
         IF ( JA .EQ. 2 .AND. JZ .NE. 1 ) THEN
            JEMISS = 1 + JZ / 2
            JRESID = JEMISS
            RNMASS = ZMASS (JRESID)
            U = 0.D+00
            DELTU = UMO - 2.D+00 * ZMASS (JEMISS)
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            IF ( DELTU .LE. 0.D+00 ) THEN
               IF ( DELTU .LT. - 2.D+00 * ANGLGB * UMO ) THEN
                  WRITE ( LUNERR, * )' *** Dres: insufficient Umo for',
     &                               ' a nucleon couple', UMO,
     &                                 2.D+00 * ZMASS (JEMISS)
               END IF
               UMO = ( UMO + DELTU ) * ( 1.D+00 + ANGLGB )
            END IF
*  |  |  |
*  |  |  +-------------------------------------------------------------*
            GO TO 2500
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
*  |  +----------------------------------------------------------------*
*  |  | Then check we are not concerning with one of the six
*  |  | standard particles
         DO 1700 J = 6, 1, -1
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            IF ( JZ .EQ. IZ (J) .AND. JA .EQ. IA (J) ) THEN
               HEVSUM = SMOM1(3) + SMOM1(5) + SMOM1(6) + SMOM1(4)
               GO TO ( 1100, 1100, 1600, 1500, 1400, 1300 ), J
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |  Proton or neutron, nothing can be done
 1100          CONTINUE
                  RETURN
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |  Alpha:
 1300          CONTINUE
                  DEUDEU = MAX ( ZERZER, U + TWOTWO * BNMASS (3)
     &                           - BNMASS (6) )
                  PROTRI = MAX ( ZERZER, U + BNMASS (4) - BNMASS (6) )
                  UEU3HE = MAX ( ZERZER, U + BNMASS (5) - BNMASS (6) )
                  QNORM  = DEUDEU + PROTRI + UEU3HE
*  |  |  |  |  If we cannot split then return
                  IF ( QNORM .LE. ZERZER ) RETURN
                  CALL GRNDM(RNDM,1)
                  V = RNDM (1)
*  |  |  |  |  +-------------------------------------------------------*
*  |  |  |  |  | Split or into two deuterons or a triton and a proton
*  |  |  |  |  | or a 3-He and a neutron: no account is made for
*  |  |  |  |  | Coulomb effects, probability is simply assumed
*  |  |  |  |  | proportional to reaction Qs
                  IF ( V .LT. DEUDEU / QNORM ) THEN
*  |  |  |  |  | Two deuterons selected
                     JEMISS = 3
                     JRESID = 3
                     RNMASS = ZMASS (3)
                     U = ZERZER
                     GO TO 2500
*  |  |  |  |  |
*  |  |  |  |  +-------------------------------------------------------*
*  |  |  |  |  | Split into a triton and a proton
                  ELSE IF ( V .LT. ( DEUDEU + PROTRI ) / QNORM ) THEN
                     JEMISS = 2
                     JRESID = 4
                     RNMASS = ZMASS (4)
                     U = ZERZER
                     GO TO 2500
*  |  |  |  |  |
*  |  |  |  |  +-------------------------------------------------------*
*  |  |  |  |  | Split into a 3-He and a neutron
                  ELSE
                     JEMISS = 1
                     JRESID = 5
                     RNMASS = ZMASS (5)
                     U = ZERZER
                     GO TO 2500
                  END IF
*  |  |  |  |  |
*  |  |  |  |  +-------------------------------------------------------*
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |  3-He:
 1400          CONTINUE
                  DEUPRO = MAX ( ZERZER, U + BNMASS (3) - BNMASS (5) )
                  PRPRNE = MAX ( ZERZER, U - BNMASS (5) )
                  QNORM  = DEUPRO + PRPRNE
*  |  |  |  |  If we cannot split then return
                  IF ( QNORM .LE. ZERZER ) RETURN
                  CALL GRNDM(RNDM,1)
                  V = RNDM (1)
*  |  |  |  |  +-------------------------------------------------------*
*  |  |  |  |  | Split or into a deuteron and a proton
*  |  |  |  |  | or into two protons and one neutron: no account is
*  |  |  |  |  | made for Coulomb effects, probability is simply assumed
*  |  |  |  |  | prportional to reaction Qs
                  IF ( V .LT. DEUPRO / QNORM ) THEN
*  |  |  |  |  | A deuteron and a proton selected
                     JEMISS = 2
                     JRESID = 3
                     RNMASS = ZMASS (3)
                     U = ZERZER
                     GO TO 2500
*  |  |  |  |  |
*  |  |  |  |  +-------------------------------------------------------*
*  |  |  |  |  | Split into 2 protons and 1 neutron: part of the exci-
*  |  |  |  |  | tation energy is conserved to allow the further
*  |  |  |  |  | splitting of the deuteron
                  ELSE
                     JEMISS = 2
                     JRESID = 0
                     FACT = ONEONE
*  |  |  |  |  |  +----------------------------------------------------*
*  |  |  |  |  |  | Loop to compute the residual excitation energy
 1450                CONTINUE
                        FACT = FACT * 0.6666666666666667D+00
*  |  |  |  |  |  | Erncm, Eepcm are the total energies of the residual
*  |  |  |  |  |  | nucleus and of the emitted particle in the CMS frame
                        U      = FACT * PRPRNE + BNMASS (3)
                        RNMASS = ZMASS (3) + U
                        ERNCM  = HLFHLF * ( UMO2 + RNMASS**2
     &                         - Z2MASS (JEMISS) ) / UMO
                        EEPCM  = UMO - ERNCM
                     IF ( EEPCM .LE. ZMASS (JEMISS) ) GO TO 1450
*  |  |  |  |  |  |
*  |  |  |  |  |  +----------------------------------------------------*
                     GO TO 2600
                  END IF
*  |  |  |  |  |
*  |  |  |  |  +-------------------------------------------------------*
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |  Triton:
 1500          CONTINUE
                  DEUNEU = MAX ( ZERZER, U + BNMASS (3) - BNMASS (4) )
                  PRNENE = MAX ( ZERZER, U - BNMASS (4) )
                  QNORM  = DEUNEU + PRNENE
*  |  |  |  |  If we cannot split then return
                  IF ( QNORM .LE. ZERZER ) RETURN
                  CALL GRNDM(RNDM,1)
                  V = RNDM (1)
*  |  |  |  |  +-------------------------------------------------------*
*  |  |  |  |  | Split or into a deuteron and a neutron
*  |  |  |  |  | or into two protons and one neutron: no account is
*  |  |  |  |  | made for Coulomb effects, probability is simply assumed
*  |  |  |  |  | proportional to reaction Qs
                  IF ( V .LT. DEUNEU / QNORM ) THEN
*  |  |  |  |  | A deuteron and a proton selected
                     JEMISS = 1
                     JRESID = 3
                     RNMASS = ZMASS (3)
                     U = ZERZER
                     GO TO 2500
*  |  |  |  |  |
*  |  |  |  |  +-------------------------------------------------------*
*  |  |  |  |  | Split into 1 proton and 2 neutrons: part of the exci-
*  |  |  |  |  | tation energy is conserved to allow the further
*  |  |  |  |  | splitting of the deuteron
                  ELSE
                     JEMISS = 1
                     JRESID = 0
                     FACT = ONEONE
*  |  |  |  |  |  +----------------------------------------------------*
*  |  |  |  |  |  | Loop to compute the residual excitation energy
 1550                CONTINUE
                        FACT = FACT * 0.6666666666666667D+00
*  |  |  |  |  |  | Erncm, Eepcm are the total energies of the residual
*  |  |  |  |  |  | nucleus and of the emitted particle in the CMS frame
                        U      = FACT * PRNENE + BNMASS (3)
                        RNMASS = ZMASS (3) + U
                        ERNCM  = HLFHLF * ( UMO2 + RNMASS**2
     &                         - Z2MASS (JEMISS) ) / UMO
                        EEPCM  = UMO - ERNCM
                     IF ( EEPCM .LE. ZMASS (JEMISS) ) GO TO 1550
*  |  |  |  |  |  |
*  |  |  |  |  |  +----------------------------------------------------*
                     GO TO 2600
                  END IF
*  |  |  |  |  |
*  |  |  |  |  +-------------------------------------------------------*
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |  Deuteron:
 1600          CONTINUE
*  |  |  |  |  +-------------------------------------------------------*
*  |  |  |  |  | Split into a proton and a neutron if it is possible
                  IF ( U .GT. BNMASS (3) ) THEN
                     JEMISS = 1
                     JRESID = 2
                     RNMASS = ZMASS (2)
                     U = ZERZER
                     GO TO 2500
*  |  |  |  |  |
*  |  |  |  |  +-------------------------------------------------------*
*  |  |  |  |  | Energy too low to split the deuteron, return
                  ELSE
                     WRITE (LUNERR,*)' **Dres: energy too low to split',
     &                               ' a deuteron! M2,M3',M2,M3
                     RETURN
                  END IF
*  |  |  |  |  |
*  |  |  |  |  +-------------------------------------------------------*
            END IF
*  |  |  |
*  |  |  +-------------------------------------------------------------*
 1700    CONTINUE
*  |  |
*  |  +----------------------------------------------------------------*
 2000    CONTINUE
         A = JA
         Z = JZ
         Q (0)  = ZERZER
         ENERG0 = FKENER (A,Z)
*  |  +----------------------------------------------------------------*
*  |  |   Note that Q(i) are not the reaction Qs but the remaining
*  |  |   energy after the reaction
         DO 2100 K = 1, 6
            JJA = JA - IA (K)
            JJZ = JZ - IZ (K)
            JJN = JJA - JJZ
            IF ( JJN .LT. 0 .OR. JJZ .LT. 0 ) THEN
               Q (K) = Q (K-1)
               GO TO 2100
            END IF
            DDJJA = JJA
            DDJJZ = JJZ
            Q (K) = MAX ( U + ENERG0 - FKENER ( DDJJA, DDJJZ )
     &            - EXMASS (K), ZERZER ) + Q (K-1)
 2100    CONTINUE
*  |  |
*  |  +----------------------------------------------------------------*
*  |  +----------------------------------------------------------------*
*  |  |  If no emission channel is open then return
         IF ( Q (6) .LE. ZERZER ) THEN
            HEVSUM = SMOM1(3) + SMOM1(5) + SMOM1(6) + SMOM1(4)
            RETURN
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         CALL GRNDM(RNDM,1)
         V = RNDM (1)
         FACT = ONEONE
*  |  +----------------------------------------------------------------*
*  |  |
         DO 2200 J = 1, 6
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            IF ( V .LT. Q (J) / Q (6) ) THEN
               JEMISS = J
               JJA    = JA - IA (JEMISS)
               JJZ    = JZ - IZ (JEMISS)
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |
               DO 2150 JJ = 1, 6
*  |  |  |  |  +-------------------------------------------------------*
*  |  |  |  |  |
                  IF ( JJA .EQ. IA (JJ) .AND. JJZ .EQ. IZ (JJ) ) THEN
                     JRESID = JJ
                     RNMASS = ZMASS (JRESID)
                     ERNCM  = HLFHLF * ( UMO2 + Z2MASS (JRESID)
     &                      - Z2MASS (JEMISS) ) / UMO
                     EEPCM  = UMO - ERNCM
                     U = ZERZER
                     GO TO 2600
                  END IF
*  |  |  |  |  |
*  |  |  |  |  +-------------------------------------------------------*
 2150          CONTINUE
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
               AJJA   = JJA
               ZJJZ   = JJZ
               RNMAS0 = AJJA * AMUMEV + FKENER ( AJJA, ZJJZ )
               GO TO 2300
            END IF
*  |  |  |
*  |  |  +-------------------------------------------------------------*
 2200    CONTINUE
*  |  |
*  |  +----------------------------------------------------------------*
         WRITE ( LUNOUT,* )' **** error in Dres, few nucleon treatment',
     &                     ' ****'
         WRITE ( LUNERR,* )' **** error in Dres, few nucleon treatment',
     &                     ' ****'
         RETURN
*  |  +----------------------------------------------------------------*
*  |  | Loop to compute the residual excitation energy
 2300    CONTINUE
            FACT = FACT * AJJA / A
            U = FACT * ( Q (JEMISS) - Q (JEMISS-1) )
*  |  | Erncm, Eepcm are the total energies of the residual
*  |  | nucleus and of the emitted particle in the CMS frame
            RNMASS = RNMAS0 + U
            ERNCM  = HLFHLF * ( UMO2 + RNMASS**2
     &             - Z2MASS (JEMISS) ) / UMO
            EEPCM  = UMO - ERNCM
         IF ( EEPCM .LE. ZMASS (JEMISS) ) THEN
            IF ( Q (JEMISS) - Q (JEMISS-1) .GE. 1.D-06 ) GO TO 2300
*  |  +--<--<--<--<--<--< Loop back
*  |  |  Actually there is no excitation energy available!
            U = ANGLGB
            RNMASS = ONEPLS * RNMAS0
            ERNCM  = ONEPLS * RNMASS
            EEPCM  = ONEPLS * ZMASS (JEMISS)
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         GO TO 2600
*  |  From here standard two bodies kinematics with Jemiss, Rnmass
*  |  (new excitation energy is U)
 2500    CONTINUE
*  |  Erncm, Eepcm are the total energies of the residual
*  |  nucleus and of the emitted particle in the CMS frame
         ERNCM = HLFHLF * ( UMO2 + RNMASS**2 - Z2MASS (JEMISS) ) / UMO
         EEPCM = UMO - ERNCM
 2600    CONTINUE
*  |  C(i) are the direction cosines of the emitted particle
*  |  (Jemiss) in the CMS frame, of course - C(i)
*  |  are the ones of the residual nucleus (Jresid if one of the
*  |  standard partcles, say the proton)
         CALL RACO (C(1),C(2),C(3))
         PCMS  = SQRT ( EEPCM**2 - Z2MASS (JEMISS) )
*  |  Now we perform the Lorentz transformation back to the original
*  |  frame (lab frame)
*  |  First the emitted particle:
         ETAX = ETACM * COSLBR (1)
         ETAY = ETACM * COSLBR (2)
         ETAZ = ETACM * COSLBR (3)
         PCMSX = PCMS * C (1)
         PCMSY = PCMS * C (2)
         PCMSZ = PCMS * C (3)
         ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
         EPS = GAMCM * EEPCM + ETAPCM - ZMASS (JEMISS)
         PHELP = ETAPCM / ( GAMCM + ONEONE ) + EEPCM
         PLBPX = PCMSX + ETAX * PHELP
         PLBPY = PCMSY + ETAY * PHELP
         PLBPZ = PCMSZ + ETAZ * PHELP
         PHELP = SQRT (PLBPX * PLBPX + PLBPY * PLBPY + PLBPZ * PLBPZ)
         COSLBP (1) = PLBPX / PHELP
         COSLBP (2) = PLBPY / PHELP
         COSLBP (3) = PLBPZ / PHELP
*  |  Then the residual nucleus ( for it c (i) --> - c (i) ):
         EREC  = GAMCM * ERNCM - ETAPCM - RNMASS
         EKRES = 1.D-03 * EREC
         PHELP = - ETAPCM / ( GAMCM + ONEONE ) + ERNCM
         PXRES = 1.D-03 * ( - PCMSX + ETAX * PHELP )
         PYRES = 1.D-03 * ( - PCMSY + ETAY * PHELP )
         PZRES = 1.D-03 * ( - PCMSZ + ETAZ * PHELP )
         P2RES = PXRES * PXRES + PYRES * PYRES + PZRES * PZRES
         PTRES = SQRT (P2RES)
         COSLBR (1) = PXRES / PTRES
         COSLBR (2) = PYRES / PTRES
         COSLBR (3) = PZRES / PTRES
*  |  Score the emitted particle
         NPART (JEMISS) = NPART (JEMISS) + 1
         SMOM1 (JEMISS) = SMOM1 (JEMISS) + EPS
         ITEMP=NPART(JEMISS)
         EPART(ITEMP,JEMISS)=EPS
         COSEVP(1,ITEMP,JEMISS)=COSLBP(1)
         COSEVP(2,ITEMP,JEMISS)=COSLBP(2)
         COSEVP(3,ITEMP,JEMISS)=COSLBP(3)
*  |  +----------------------------------------------------------------*
*  |  |  Check if the residual nucleus is one of the emitted particles
         IF ( JRESID .GT. 0 ) THEN
            J = JRESID
            GO TO 6102
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         JA = JA - IA (JEMISS)
         JZ = JZ - IZ (JEMISS)
* Umo is the invariant mass of the system!!
         UMO  = RNMASS
         UMO2 = UMO * UMO
         ELBTOT  = RNMASS + EREC
         GAMCM   = ELBTOT / RNMASS
         ETACM   = 1.D+03 * PTRES / RNMASS
         GO TO 1000
      END IF
*  |
*  +-------------------------------------------------------------------*
* Come to 23 at the beginning and after the end of a "normal"
* evaporation cycle
   23 CONTINUE
      A = JA
      Z = JZ
      IF(JA-8)8486,8488,8486
 8488 CONTINUE
      IF(JZ-4)8486,1224,8486
 8486 CONTINUE
      DO 2 K=1,6
         IF((A-FLA(K)).LE.(Z-FLZ(K))) THEN
            Q(K)=99999.D0
            GO TO 2
         END IF
         IF(A-TWOTWO*FLA(K))2,727,727
  727    CONTINUE
         IF(Z-TWOTWO*FLZ(K))2,728,728
  728    CONTINUE
         Q(K) = QNRG(A-FLA(K),Z-FLZ(K),A,Z) + EXMASS(K)
C 728    Q(K)=FKENER(A-FLA(K),Z-FLZ(K))-FKENER(A,Z)+EXMASS(K)
    2 CONTINUE
      FLKCOU(1)=ZERZER
      FLKCOU(2)=DOST(1,Z-FLZ(2))
      FLKCOU(3)=FLKCOU(2)+.06D+00
      FLKCOU(4)=FLKCOU(2)+.12D+00
      FLKCOU(6)=DOST(2,Z-FLZ(6))
      FLKCOU(5)=FLKCOU(6)-.06D+00
      CCOUL(1)=ONEONE
      CCOU2=DOST(3,Z-FLZ(2))
      CCOUL(2)=CCOU2+ONEONE
      CCOUL(3)=CCOU2*1.5D0+THRTHR
      CCOUL(4)=CCOU2+THRTHR
      CCOUL(6)=DOST(4,Z-FLZ(6))*TWOTWO+TWOTWO
      CCOUL(5)=TWOTWO*CCOUL(6)-ONEONE
      SIGMA=ZERZER
*  Initialize the flag which checks for open particle decay with
*  zero excitation and pairing --> for particle unstable residual
*  nuclei
      LOPPAR = .FALSE.
      DO 33 J=1,6
         IF(A-TWOTWO*FLA(J))5,725,725
  725    CONTINUE
         IF(Z-TWOTWO*FLZ(J))5,726,726
  726    CONTINUE
         MM=JA-IA(J)
         ZZ=Z-FLZ(J)
         AA=A-FLA(J)
         IF(AA.LE.ZZ)GO TO 5
*  Energy threshold for the emission of the jth-particle
         THRESH(J)=Q(J)+.88235D+00*FLKCOU(J)*FLZ(J)*
     &   ZZ/(RMASS(MM)+RHO(J))
         LOPPAR = LOPPAR .OR. THRESH (J) .GT. ZERZER
         IAA = NINT(AA)
         IZZ = NINT(ZZ)
         NN  = IAA - IZZ
*  The residual nucleus excitation energy ranges from 0 up
*  to U - Q (J)
         ILVMOD = IB0
         UMXRES = U - THRESH (J)
*  This is the a lower case of the level density
         SMALLA (J) = GETA ( UMXRES, IZZ, NN, ILVMOD, ISDUM, ASMMAX,
     &                       ASMMIN )
         CALL EEXLVL (IAA,IZZ,EEX1ST,EEX2ND,CORR)
         EEX1ST = 1.D+03 * EEX1ST
         EEX2ND = 1.D+03 * EEX2ND
         CORR   = 1.D+03 * MAX ( CORR, ZERZER )
         IF ( NN .EQ. 4 .AND. IZZ .EQ. 4 ) THEN
            IF ( U - THRESH (J) - 6.1D+00 .GT. ZERZER ) THEN
               CORR = SIXSIX
            ELSE
               TMPVAR = U-THRESH(J)-0.1D+00
               CORR = MAX ( ZERZER, TMPVAR )
            END IF
         END IF
         IF (NINT(FKEY).EQ.1) CORR=ZERZER
         CORRRR(J)=CORR
* Standard calculation:
         ARG=U-THRESH(J)-CORR
         IF(ARG)5,4,4
    5    CONTINUE
         R(J)=ZERZER
         S(J)=ZERZER
         SOS(J)=ZERZER
         GO TO 33
    4    CONTINUE
         S(J)=SQRT (SMALLA(J)*ARG)*TWOTWO
         SOS(J)=TENTEN*S(J)
   33 CONTINUE
      N1=1
      SES = MAX (S(1),S(2),S(3),S(4),S(5),S(6))
      DO 1131 J=1,6
         JS  = SOS(J) + ONEONE
         FJS = JS
         STRUN(J) = FJS - ONEONE
         IF ( S(J) .GT. ZERZER ) THEN
            MM  = JA-IA(J)
            EXPSAS=MIN(EXPMAX,MAX(EXPMIN,S(J)-SES))
            SAS = EXP (EXPSAS)
            EXPSUS=MIN(EXPMAX,MAX(EXPMIN,-S(J)))
            SUS = EXP (EXPSUS)
            EYE1(J) = ( TWOTWO * S(J)**2 -SIXSIX * S(J)
     &              + SIXSIX + SUS * ( S(J)**2 - SIXSIX ) )
     &              / ( EIGEIG * SMALLA(J)**2 )
            IF ( J .EQ. 1 ) THEN
               EYE0(J) = ( S(J) - ONEONE + SUS ) / ( TWOTWO*SMALLA(J) )
* Standard calculation
               R   (J) = RMASS(MM)**2 * ALPH(MM) * ( EYE1(J) + BET(MM)
     &                 * EYE0(J) ) * SAS
            ELSE
               R   (J) = CCOUL(J) * RMASS(MM)**2 * EYE1(J) * SAS
            END IF
            R (J) = MAX ( ZERZER, R (J) )
            SIGMA = SIGMA + R (J)
         END IF
 1131 CONTINUE
      NCOUNT = 0
 6202 CONTINUE
      IF(SIGMA)9,9,10
 9    CONTINUE
      DO 6100 J = 1,6
         IF(JA-IA(J))6100,6101,6100
 6101    CONTINUE
         IF(JZ-IZ(J))6100,6102,6100
 6100 CONTINUE
      GO TO 6099
 6102 CONTINUE
      IF ( U .GT. ANGLGB ) GO TO 1000
      JEMISS = J
C*****STORE,RESIDUAL NUC IS OF EMITTED PARTICLE TYPE
* If we are here this means that the residual nucleus is equal to
* one of the six emitted particle (the j-th one). So give to it
* all the energy, score it and return with 0 recoil and excitation
* energy for the residual nucleus
      EPS = EREC
      NPART(JEMISS) = NPART(JEMISS)+1
      ITEMP=NPART(JEMISS)
      NRNEEP = NRNEEP + 1
      SMOM1(JEMISS) = SMOM1(JEMISS) + EPS
      ITEMP=NPART(JEMISS)
      EPART(ITEMP,JEMISS)=EPS
      COSEVP(1,ITEMP,JEMISS) = COSLBR(1)
      COSEVP(2,ITEMP,JEMISS) = COSLBR(2)
      COSEVP(3,ITEMP,JEMISS) = COSLBR(3)
      GO TO 8002
*
*-->-->-->-->--> go to return
 6099 CONTINUE
      IF(JA-8)72,51,72
   51 CONTINUE
      IF(JZ-4)72,1224,72
* Come here for a "normal" step
   10 CONTINUE
      LOPPAR = .FALSE.
      CALL GRNDM(RNDM,1)
      URAN=RNDM(1)*SIGMA
      SUM = ZERZER
      DO 16 J=1,6
         K   = J
         SUM = R(J)+SUM
         IF ( SUM - URAN .GT. 0.D+00 ) GO TO 17
   16 CONTINUE
   17 CONTINUE
      JEMISS=K
      NPART(JEMISS) = NPART (JEMISS) + 1
      JS = SOS (JEMISS) + ONEONE
      IF(JS-1000)4344,4345,4345
 4345 CONTINUE
      RATIO2=(S(JEMISS)**3-6.D0*S(JEMISS)**2+15.D0*
     &S(JEMISS)-15.D0)/((2.D0*S(JEMISS)**2-6.D0*S(JEMISS)+6.D0)*SMALLA
     &(JEMISS))
      GO TO 4347
 4344 CONTINUE
      RATIO2=(P2(JS)+(P2(JS+1)-P2(JS))*
     &(SOS(JEMISS)-STRUN(JEMISS)))/SMALLA(JEMISS)
 4347 CONTINUE
      EPSAV=RATIO2*TWOTWO
*  +-------------------------------------------------------------------*
*  |  Neutron channel selected:
      IF (JEMISS .EQ. 1) THEN
         MM=JA-IA(J)
         EPSAV=(EPSAV+BET(MM))/(ONEONE+BET(MM)*EYE0(JEMISS)
     &        /EYE1(JEMISS))
*  |  +----------------------------------------------------------------*
*  |  |  Compute the fission width relative to the neutron one:
*  |  |  this part is taken from subroutine EMIT of LAHET
         IF ( IFISS .GT. 0 .AND. JZ .GE. 71 .AND. .NOT. FISINH ) THEN
*  |  |  +-------------------------------------------------------------*
*  |  |  |  Compute the correction factor for the fission width:
            IF ( JZ .GT. 88 ) THEN
               AGOES = ONEONE
*  |  |  |
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            ELSE
               AGOES = MAX ( ONEONE, ( U-SEVSEV ) / ( EPSAV+SEVSEV ) )
            END IF
*  |  |  |
*  |  |  +-------------------------------------------------------------*
*  |  |  Finally this is the relative fission width:
*  |  |  This is : Probfs = 1 / ( 1 + G_n / G_f )
            PROBFS = FPROB ( Z, A, U ) / AGOES
*  |  |  +-------------------------------------------------------------*
*  |  |  |  Check if it will be fission:
            CALL GRNDM(RNDM,1)
            IF ( RNDM (1) .LT. PROBFS ) THEN
               FISINH = .TRUE.
               KFISS  = 1
*  |  |  |  Undo the counting of the neutron evaporation
               NPART (JEMISS) = NPART (JEMISS) -1
               GO TO 9260
            END IF
*  |  |  |
*  |  |  +-------------------------------------------------------------*
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
      END IF
*  |
*  +-------------------------------------------------------------------*
  20  CONTINUE
      CALL GRNDM(RNDM,2)
      E1=-HLFHLF*LOG(RNDM(1))
      E2=-HLFHLF*LOG(RNDM(2))
* Eps should be the total kinetic energy in the CMS frame
* Standard calculation:
      EPS=(E1+E2)*EPSAV+THRESH(JEMISS)-Q(JEMISS)
      AR = A - IA(JEMISS)
      ZR = Z - IZ(JEMISS)
* The CMS energy is updated
      IMASS = NINT (AR)
      IF ( IMASS .EQ. 8 .AND. NINT (ZR) .EQ. 4 ) THEN
         UNEW = U - EPS - Q(JEMISS)
         UMAX = U - THRESH(JEMISS)
         IF ( UNEW .GT. 6.D+00 ) THEN
            UMIN = 6.D+00
         ELSE IF ( UNEW .GT. 4.47D+00 .AND. UMAX .GT. 6.D+00 ) THEN
            UMIN = 4.47D+00
            UNEW = 6.D+00
         ELSE IF ( UNEW .GT. 1.47D+00 .AND. UMAX .GT. 2.94D+00 ) THEN
            UMIN = 1.47D+00
            UNEW = 2.94D+00
         ELSE
            UMIN = -0.1D+00
            UNEW = ANGLGB * 0.1D+00
         END IF
      ELSE IF ( IMASS .LE. 4 ) THEN
         IPRO = NINT ( ZR )
         INEU = IMASS - IPRO
         IF ( IMASS .EQ. 1 ) THEN
*  Be sure that residual neutrons or protons are not left excited
            UMIN = 0.D+00
            UNEW = 0.D+00
            EPS  = U - Q(JEMISS)
         ELSE IF ( IPRO .EQ. 0 .OR. INEU .EQ. 0 ) THEN
*  Ipro protons or ineu neutrons arrived here!
            UMIN = CORRRR(JEMISS)
            UNEW = U - EPS - Q(JEMISS)
         ELSE IF ( IMASS .LE. 2 ) THEN
*  Be sure that residual deuterons are not left excited!
            UMIN = 0.D+00
            UNEW = 0.D+00
            EPS  = U - Q(JEMISS)
         ELSE IF ( ABS ( INEU - IPRO ) .LE. 1 ) THEN
*  For the moment also residual 3-H, 3-He and 4-He are not left
*  excited !
            UMIN = 0.D+00
            UNEW = 0.D+00
            EPS  = U - Q(JEMISS)
         ELSE
            UMIN = CORRRR(JEMISS)
            UNEW = U - EPS - Q(JEMISS)
         END IF
      ELSE
         UMIN = CORRRR(JEMISS)
         UNEW = U - EPS - Q(JEMISS)
      END IF
* Standard calculation
      IF(UNEW-UMIN)6200,6220,6220
 6220 CONTINUE
*or   RNMASS = AR * AMUMEV + FKENER (AR,ZR)
      RNMASS = AR * AMUMEV + FKENER (AR,ZR) + UNEW
      UMIN2  = ( RNMASS + ZMASS (JEMISS) )**2
      IF ( UMIN2 .GE. UMO2 ) THEN
         GO TO 6200
      END IF
      U = UNEW
* C(i) are the direction cosines of the evaporated particle in the CMS
* frame, of course - C(i) are the ones of the residual nucleus
      CALL RACO(C(1),C(2),C(3))
* Erncm, Eepcm are the total energies of the residual nucleus and
* of the evaporated particle in the CMS frame
      ERNCM = 0.5D+00 * ( UMO2 + RNMASS**2 - Z2MASS (JEMISS) ) / UMO
      EEPCM = UMO - ERNCM
      PCMS  = SQRT ( EEPCM**2 - Z2MASS (JEMISS) )
* Now we perform the Lorentz transformation back to the original
* frame (lab frame)
* First the evaporated particle:
      ETAX = ETACM * COSLBR (1)
      ETAY = ETACM * COSLBR (2)
      ETAZ = ETACM * COSLBR (3)
      PCMSX = PCMS * C (1)
      PCMSY = PCMS * C (2)
      PCMSZ = PCMS * C (3)
      ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
      EPS = GAMCM * EEPCM + ETAPCM - ZMASS (JEMISS)
      PHELP = ETAPCM / (GAMCM + 1.D0) + EEPCM
      PLBPX = PCMSX + ETAX * PHELP
      PLBPY = PCMSY + ETAY * PHELP
      PLBPZ = PCMSZ + ETAZ * PHELP
      PHELP = SQRT (PLBPX * PLBPX + PLBPY * PLBPY + PLBPZ * PLBPZ)
      COSLBP (1) = PLBPX / PHELP
      COSLBP (2) = PLBPY / PHELP
      COSLBP (3) = PLBPZ / PHELP
* Then the residual nucleus ( for it c (i) --> - c (i) ):
      EREC  = GAMCM * ERNCM - ETAPCM - RNMASS
      EKRES = 1.D-03 * EREC
      PHELP = - ETAPCM / (GAMCM + 1.D0) + ERNCM
      PXRES = 1.D-03 * ( - PCMSX + ETAX * PHELP )
      PYRES = 1.D-03 * ( - PCMSY + ETAY * PHELP )
      PZRES = 1.D-03 * ( - PCMSZ + ETAZ * PHELP )
      P2RES = PXRES * PXRES + PYRES * PYRES + PZRES * PZRES
      PTRES = SQRT (P2RES)
      COSLBR (1) = PXRES / PTRES
      COSLBR (2) = PYRES / PTRES
      COSLBR (3) = PZRES / PTRES
* Check energy and momentum conservation !!
      IF (EREC .LE. 0.D+00) THEN
         PTRES = 0.D+00
         EREC  = 0.D+00
      END IF
* Umo is the invariant mass of the system!!
      UMO  = RNMASS
      UMO2 = UMO * UMO
      ELBTOT  = RNMASS + EREC
      GAMCM   = ELBTOT / RNMASS
      ETACM   = 1.D+03 * PTRES / RNMASS
      GO TO 76
 6200 CONTINUE
      NCOUNT = NCOUNT + 1
      IF ( NCOUNT .LE. 10 ) GO TO 20
      SIGMA = SIGMA - R(JEMISS)
* if we are here we have sampled for > 10 times a negative energy Unew
      NPART(JEMISS)=NPART(JEMISS)-1
      R(JEMISS) = 0.D0
      NCOUNT = 0
      GO TO 6202
   76 CONTINUE
      JAT=JA-IA(JEMISS)
      JZT=JZ-IZ(JEMISS)
      IF(JAT-JZT)9,9,77
   77 CONTINUE
      JA=JAT
      JZ=JZT
C*****STORE,END OF NORMAL CYCLE
      SMOM1(JEMISS)=SMOM1(JEMISS)+EPS
      ITEMP=NPART(JEMISS)
      EPART(ITEMP,JEMISS)=EPS
      COSEVP(1,ITEMP,JEMISS)=COSLBP(1)
      COSEVP(2,ITEMP,JEMISS)=COSLBP(2)
      COSEVP(3,ITEMP,JEMISS)=COSLBP(3)
* The following card switch to the rough splitting treatment
      IF (JA .LE. 2) GO TO 1000
      IF(JA-8)23,1223,23
 1223 CONTINUE
      IF(JZ-4)23,1224,23
* If we are here the residual nucleus is a 8-Be one, break it into
* two alphas with all the available energy (U plus the Q of the breakup)
* , score them and return with 0 recoil and excitation energy
 1224 CONTINUE
      LOPPAR = .FALSE.
      IF(U)1228,1229,1229
 1228 CONTINUE
      EPS=0.D0
      GO TO 1230
 1229 CONTINUE
 1230 CONTINUE
      NBE8=NBE8+1
* C(i) are the direction cosines of the first alpha in the CMS
* frame, of course - C(i) are the ones of the other
      CALL RACO(C(1),C(2),C(3))
* Ecms is the total energy of the alphas in the CMS frame
      ECMS  = 0.5D+00 * UMO
      PCMS  = SQRT ( ECMS**2 - Z2MASS (6) )
* Now we perform the Lorentz transformation back to the original
* frame (lab frame)
* First alpha:
      ETAX = ETACM * COSLBR (1)
      ETAY = ETACM * COSLBR (2)
      ETAZ = ETACM * COSLBR (3)
      PCMSX = PCMS * C (1)
      PCMSY = PCMS * C (2)
      PCMSZ = PCMS * C (3)
      ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
      EPS = GAMCM * ECMS + ETAPCM - ZMASS (6)
      PHELP = ETAPCM / (GAMCM + 1.D0) + ECMS
      PLBPX = PCMSX + ETAX * PHELP
      PLBPY = PCMSY + ETAY * PHELP
      PLBPZ = PCMSZ + ETAZ * PHELP
      PHELP = SQRT (PLBPX * PLBPX + PLBPY * PLBPY + PLBPZ * PLBPZ)
* Store the first alpha!!
      SMOM1(6) = SMOM1(6) + EPS
      NPART(6) = NPART(6) + 1
      ITEMP = NPART(6)
      EPART (ITEMP,6) = EPS
      COSEVP (1,ITEMP,6) = PLBPX / PHELP
      COSEVP (2,ITEMP,6) = PLBPY / PHELP
      COSEVP (3,ITEMP,6) = PLBPZ / PHELP
* Then the second alpha ( for it c (i) --> - c (i) ):
      EPS = GAMCM * ECMS - ETAPCM - ZMASS (6)
      PHELP = - ETAPCM / (GAMCM + 1.D0) + ECMS
      PLBPX = - PCMSX + ETAX * PHELP
      PLBPY = - PCMSY + ETAY * PHELP
      PLBPZ = - PCMSZ + ETAZ * PHELP
      PHELP = SQRT (PLBPX * PLBPX + PLBPY * PLBPY + PLBPZ * PLBPZ)
* Store the second alpha !!
      SMOM1(6) = SMOM1(6) + EPS
      NPART(6) = NPART(6) + 1
      ITEMP = NPART(6)
      EPART (ITEMP,6) = EPS
      COSEVP (1,ITEMP,6) = PLBPX / PHELP
      COSEVP (2,ITEMP,6) = PLBPY / PHELP
      COSEVP (3,ITEMP,6) = PLBPZ / PHELP
 8002 CONTINUE
      LOPPAR = .FALSE.
      EREC   = ZERZER
      U      = ZERZER
      EKRES  = ZERZER
      PTRES  = ZERZER
   72 CONTINUE
      HEVSUM=SMOM1(3)+SMOM1(5)+SMOM1(6)+SMOM1(4)
      RETURN
*.......................................................................
*///// RAL FISSION ROUTINE /////
 9260 CONTINUE
*  +-------------------------------------------------------------------*
*  |  Record the direction cosines of the fissioning nucleus
      DO 9270 I=1,3
         COSLF0 (I) = COSLBR (I)
 9270 CONTINUE
*  |
*  +-------------------------------------------------------------------*
      CALL FISFRA ( JA, JZ, U, EREC, UMO, GAMCM, ETACM )
*  +-------------------------------------------------------------------*
*  |  Check for fission failures!!
      IF ( .NOT. FISINH ) THEN
         PENBAR = .FALSE.
         GO TO 23
      END IF
*  |
*  +-------------------------------------------------------------------*
*  Do not pick up the fission fragments, rather go back to Evevap
      HEVSUM=SMOM1(3)+SMOM1(5)+SMOM1(6)+SMOM1(4)
      RETURN
      END
